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Abstract 

We present a computer simulation study, via lattice Boltzmann simulations, of a microscopic model for cy- 
toplasmic streaming in algal cells such as those of Chara corallina. We modelled myosin motors tracking along 
actin lanes as spheres undergoing directed motion along fixed lines. The sphere dimension takes into account 
the fact that motors drag vesicles or other organelles, and, unlike previous work, we model the boundary close 
to which the motors move as walls with a finite slip layer. By using realistic parameter values for actin lane and 
myosin density, as well as for endoplasmic and vacuole viscosity and the slip layer close to the wall, we find 
that this simplified view, which does not rely on any coupling between motors, cytoplasm and vacuole other than 
that provided by viscous Stokes flow, is enough to account for the observed magnitude of streaming velocities in 
intracellular fluid in living plant cells. 

Keywords: cytoplasmic streaming, lattice Boltzmann, simulation, wall slip 



Introduction 



Within eukaryotic cells, intracellular flows are often unremarkable: typical Reynolds numbers are low, and cells 
are tiny enough to quench most fluid motion within them. An important exception is provided by cytoplasmic 
streaming, which is the name given to the ~ 80 — 100 fim/s directed flow of cytosol and organelles around large 
plant and fungal cells Hl-d] 0- Unlike typical mammalian cells, which are tens of microns in size, the internode 
stalks of the Chara corallina algae, for instance, are single cells which are about a millimeter in diameter and 
several centimeters long. Giant single cells like these need to overcome a highly non-trivial transport problem in 
order to deliver nutrients and other chemicals throughout their interior. Given such length scales, thermal diffusion 
does not provide a viable option to move things around quickly enough: for instance, given typical intracellular 
diffusion coefficients, it would take months even for highly mobile ions to traverse the cell length. Cytoplasmic 
streaming gives a very efficient alternative, as it sets the whole fluid within the cell in motion, thereby readily 
advecting proteins and nutrients which happen to be there - in other words while the Reynolds number may be 
infinitesimally small, the Peclet number is instead rather large flQ]. 

Understanding the basic biophysical mechanism underpinning cytoplasmic streaming has been the subject of 
an ongoing debate at the interface between fluid dynamics and cell biology. To introduce the competing scenarios 
proposed in the literature, it is useful to recall the geometry of algal cells such as Chara corallina. The basic 
components are the "internodes" or internode stalks, which as mentioned above ai^e extra-long single cells. Each of 
these internodes has an approximately cylindrical symmetry, and may be thought of as consisting of two concentric 
layers: the external one is a micron-thick cellular layer, filled with cytosol (endoplasm) and separated from the 
large interior, known as vacuole, by a fluid membrane. The vacuole is made up of an aqueous solvent, and contains 
no internal structure. The inner surface of the external wall of the internode cell is instead covered by actin bundles, 
which act as conveyor belts along which myosin motors walk. The actin bundles, or lanes, follow helical paths on 
the cylindrical walls of Chara, and they are organised into alternating bands having opposite directions of velocity. 
The bands are separated by a so-called indifferent zone. The hypothesised mechanism for cytoplasmic streaming 
is that the directed motion of myosin motors on the actin bundle tracks can somehow entrain the cytoplasmic 
fluid, and that this in turn sets the vacuole into motion by transferring momentum through the fluid membrane 
separating it from the cytoplasm. According to this view, the myosin traffic provides an effective shear boundary 
condition on a Stokes equation for the velocity inside the vacuole. The fact that the tracks along which motors 
move are helicoidal leads to an interesting fluid dynamics problem. Its solution, presented in includes a 
secondary flow along the planes perpendicular to the cylindrical axis, which helps nutrient and chemical mixing 
along the flow gradient direction as well - mixing is notoriously a hard requirement to fulfil at zero Reynolds 
number. More recently, direct NMR velocimetry experiments of flow inside the vacuole of single intemodal cells 
have quantitatively validated this analytical solution 

Earlier experiments on cytoplasmic streaming based on laser light scattering experiments on M^eZZaJslQ] and 
Chara ifiol] cells showed that the fluid velocity distribution is quite nan^ow. While Sattelle & Buchan um speak 
of plug-like motion for both Chara and Nitella, Mustacich & Ware |@] underscore the fact that the peak in the 
spectrum of scattered light is not sharp enough for the streaming flow of a plug but indicates that most particle 
velocities are within 10% of the most likely velocity. Experiments agree that it is unlikely that any particles move 
at substantially larger velocities than those caiTied by cytoplasmic streaming - in other words, the fluid velocity 
has to be very close to that of the vesicles dragged along by the motors and driving the flow. 

The idea that myosin motion can provide an effectively smooth boundary condition for the fluid dynamics 
inside the vacuole by necessity requires a strong hydrodynamic coupling between actin-myosin motion on the cell 
surface and fluid flow in the vacuole. Finding exactly what structure can provide this coupling is however a non- 
trivial task. In an interesting contribution lllin . Nothnagel and Webb analysed the hydrodynamic feasibility of three 
models for momentum transfer from myosin to endoplasm and vacuole. Their calculation shows that individual 



'Cytoplasmic streaming may also be observed in non-plant cells, for instance in the nematode Caenorhabditis elegans Ql- 



myosin molecules running on the actin tracks are by themselves ineffective in setting the cytosol or the vacuole 
into motion. A second model considered attachment of myosin to organelles and vesicles in the endoplasm: while 
this much improves the viscous coupling, the calculations in Ref. jllll still suggested that if this mechanism was 
at the origin of cytoplasmic streaming it would require very packed traffic of motors along the actin lanes, which 
is unlikely in the real system. The final model, favoured by the semianalytical estimates in Ref. jllh . envisages 
the existence of an elastic network, or gel, incorporating the moving motors and extending into the endoplasm: in 
this framework the movement of the motors pulls the network forward in a plug-like fashion and the coupling to 
the vacuole is readily achieved. 

Subsequent experiments (video-enhanced light microscopy and fast-freezing electron microscopy) further sug- 
gested that the myosin motors may actually be attached to the endoplasmic reticulum which then performs a sliding 
motion over the actin cables il2l] . In some later papers this picture has been adopted 11311 . while others are less 
definite and speak only of cargoes in general il4il or state that the organelles associated with myosins have not 
been identified L2il- 

Here we revisit the view that the myosin motors need not be directly associated with a network structure but 
that the high viscosity of the cytoplasm together with a thin slip layer are sufficient to make it move with a very flat 
velocity gradient at roughly the same velocity as the active particles. We perform lattice Boltzmann simulations of 
micron-sized spherical vesicles moving close to a planar wall (dragged by the motors) and show that, depending on 
the slip allowed at the wall, this simpler picture can explain the uniform streaming profile also for realistic densities 
of spherical particles. The effect of the endoplasmic reticulum and all other cell contents is incorporated into the 
high viscosity of the cytoplasm, but our results suggest that it is not necessary to invoke a physical tethering of the 
motors to any network in order to explain the known phenomenology of cytoplasm streaming. 



Models and methods 

Lattice Boltzmann simulations of colloidal spheres in a fluid 

Rather than address the cylindrical geometry suggested by the Chara system, we address here a simpler, planar 
model (Fig. [Hi. This allows us to address the issues of principle outlined above while avoiding the numerical 
difficulties of a curved geometry. Moreover, as detailed below, the physics of interest to us concerns the cytosol 
and only a thin region of the vacuole containing the endoplasm-vacuole interface: this region is almost planar. 
We consider a slab of material made up of two layers with different viscosities, r// and rju > Vi' corresponding 
to respectively the vacuole and the endoplasm, modelled here as a Newtonian fluid with high viscosity, as in the 



semianalytical treatment in Ref. Bill] . In each of the layers, the fluid obeys the Navier-Stokes equation, 



p {dt + Updp) Ua = -daP + dp {daUp + riidpUa) (1) 

where p denotes fluid density, u fluid velocity, p pressure, rji is the fluid viscosity, and i = I,II labels the layer 
under consideration. In our formalism Greek letters denote Cartesian components and summation over repeated 
indices is implied. In order to solve Eq.[Tl we employ our Ludwig code, which is based on a lattice Boltzmann (LB) 



algorithm (see Ref. 111511 for details). Briefly, the LB method is based on a set of mesoscopic velocity distribution 
functions, /j(x, t), which depend on some lattice velocity vectors e^. In our 3D case, we consider 19 velocity 
vectors (see Appendix), this is known as a D3Q19 lattice. The density and velocity fields of the fluid may be 
recovered as moments of such distribution functions, as follows, 

p(x,t) = J]/i(x,t) (2) 

i 

/)(x,t)iiQ,(x,t) = y2fi{x,t)ei^a. 



In our LB algorithm, the distribution functions evolve according to the following discretised, or lattice, Boltzmann 
equation: 

/,(x + e,At, t + At) = /,(x, t) + (3) 

T 

where At is the time step (normally unity in LB), r is a relaxation parameter, related to the fluid viscosity t] 
via 1] = |(r — Ai/2), and f^'^ is a set of equilibrium distribution functions. Note that r differs in the vacuole 
and cytoplasm as the viscosities differ. The macroscopic equations of motion obeyed by p and Ua may be found 
by taking moments of Eq. [3] and expanding for small At - a. formal procedure known as the Chapman-Enskog 
expansion. The LB method defines some dynamical rules for the distribution functions which correspond to a 
continuity equation for the fluid density and Eq. \T\ for its velocity field. It can be shown that this is the case 
provided the following constraints on the moments of the equilibrium distribution are satisfied: 

J]/r(x,t) = Kx,t) (4) 

i 

i 

^fi'^{^,t)ei,aei,i3 = p(x,t)na(x,t)n/3(x,t). 

i 

At this point, it should be noted that any set of equilibrium distribution functions satisfying these constraints leads 
to a viable LB algorithm to solve the fluid dynamics problem we are interested in. A popular choice in the literature 
is to expand the equilibrium distribution functions as a power series in the velocity vectors, with the coefficients 
being determined, in general in a non-unique way, through the constraints. Note that as we want to model a two- 
layered fluid, made up by a more viscous endoplasm and a less viscous vacuole, we considered a generalisation of 
the LB equation given above, to allow for a spatially dependent value of the relaxation constant r. We found this 
simple procedure to describe well our two-layered system - we anticipate this may not hold in more complicated 
systems where the interface is not planai\ 

The LB method, in the formulation which we adopt, has been developed to allow consistent coupling between 
the fluid dynamics of a Newtonian fluid and the dynamics of spherical colloidal particles embedded in it. The 
coupling is through the well-known method of bounce-back on links, which accurately enforces no-slip boundary 
conditions for the velocity field on the surface of each particle. Instead of repeating the discussion of these 
colloidal boundary conditions here, we refer the interested reader to 1 15, 13] for details. In our case, the solid 



particles represent the motors with their associated cargoes (vesicles and organelles). They are positioned at a 
fixed distance from the bottom wall, and are subjected to a strong harmonic potential which virtually confines 
their motion to fixed ID lines, which represent the actin bundles. 

A noteworthy advantage of the LB methodology we use is that, as it fully handles the boundary conditions 
for the velocity fields both on the particle surface and on the walls, both near and far field effects are correctly 
included. With respect to e.g. Stokesian dynamics, which would be another perfectly valid choice for the problem 
at hand, LB requires one to introduce both an inertial term and a finite compressibility: both these terms need to 
be kept small so that they do not effectively contribute to the physics of our system. On the other hand, LB allows 
one to reach much larger systems than Stokesian dynamics, at least for the nonuniform geometry addressed here. 
Finally, we note that an analytical approximation based on modelling the vesicles as Stokeslets near a wall would 
not incorporate near-field interactions which are important in our simulations, as particles may come into close 
contact. 



Boundary conditions for slip and no-slip walls 



In this Section we discuss in some detail how to implement slip boundary conditions at the bottom wall (see 
geometry in Fig. 1), which is important for our work. Here we outline the main equations which we use in our 




Figure 1: (Online version in colour.) Cartoon illustration of spherical vesicles moving through the cytoplasm 
(dragged by molecular motors not shown). The cartoon is roughly to scale with vesicle diameters 1 //m, height 
of the cytoplasm 10 //m and thickness of actin bundles 100 nm. The light blue tonoplast membrane separates 
cytoplasm and vacuole which extends beyond the figure. In the simulations, the motors and the membrane are not 
modelled explicitly. 



algorithm, whereas a full derivation is given in the Appendix. 

In order to model slip at a solid wall in LB simulations, a partial slip parameter p can be introduced which 
intuitively corresponds to the "fraction of slip" at the wall. More accurately, p is the fraction of each of the 
distribution functions which is "bounced forward" at the wall, whereas 1 — p is the fraction which is "bounced 
back" (see Appendix). The limit of p — > corresponds to the commonly employed no slip boundary condition. 



The amount of wall slip may also be related to the Knudsen number of the fluid 11171] . In our scheme, we use it as 
a parameter to be determined, in practice, from the knowledge of the shp length, u(0)/n'(0), which we assume 
comes from either experiment or theory. 

Our choice of a slip wall for the bottom boundary is motivated by the need to model the fluid dynamics in 
the endoplasm close to the dragged cargo. We assume that below the cargo lies a thin aqueous low-viscosity 
layer, as the viscosity in the endoplasm derives from the presence of a labile reticular network with mesh size 
smaller than the distance between sphere and wall, or from macromolecular crowding in the endoplasm. Either 
way, the network or crowding agents, which ai^e coarse-grained into a high-viscosity fluid for the bulk endoplasm, 
ai"e depleted in the vicinity of the wall, and this accounts for the thin layer of lower viscosity. If the behaviour 
of the thin layer itself is not of interest, all it does is to introduce effective slip for the main system which can be 
simulated using the slip parameter p. 

To make progress, we start by writing down the velocity profile of the fluid as a power series in the distance 
from the bottom wall z, 

^ ^' < z < e 

e<z<R 




The slip length u{0)/u'{0) can be determined analytically for e <C i? where R is the total system size and e the 
thickness of the low-viscosity layer. By enforcing no slip at the bottom wall, together with continuity of the stress 
and velocity fields at 2 = e, and subsequently taking the limit e ^ (see Appendix), we obtain the following 
form for the slip length, 

^(0) _Jm {\ 



u'{0) V VI 



to first order in e, where r]i is the (lower) viscosity of the thin layer close to the wall and rju the (higher) viscosity 
in the bulk of the system. Here we assume that the viscosity of the thin layer is the same as that of the vacuole 
as both are basically aqueous fluids. This is however only a conceptual simplification and not necessary for the 
derivation of the boundary condition. 

The slip length u{0)/u'{0) and slip pai^ameter p may be shown to be related by the following formula (derived 
in the Appendix): 

u{0) _ (2r - 1) P _ „ P ... 

The last two equations also lead to the following explicit formula for p, 

P = ^^^^^7^^^ (7) 

^ViVii + e [1111 - Vi) 

As demonstrated in the Appendix, numerical simulations of simple test cases performed with a given p lead indeed 
to a flow profile that agrees with our analytical calculation. In the formulas above, all quantities are given in terms 
of lattice, or simulation, units. To convert these to physical units, e has to be multiplied by the lattice spacing Ax, 
while r] should be multiplied by the fluid density p and by Ax^, and divided by the time step At. In lattice, or 
simulation units, p = Ax = At = 1. 

In physical units we choose Ax = 0.4 pm, At = 0.08 ms, and p ~ 2.5 x 10^ kg/m^. With these values, 
our vesicle size and vacuole viscosity map to 0.5pm and IcP respectively, as in the experiments. Furthermore a 
velocity of 0.01 in lattice units, typically used in our simulations, maps to a velocity of 50 pm/s. However our 
chosen density is much higher than that of waterjl 

Choice of parameters 

Simulations were run for colloidal particles of radius cq = 1-25 lattice units (LU), whereas the size of the cy- 



toplasm was varied within a physically meaningful regime. In keeping with the parameters of ref. Ill 111 vesicles 
are assumed to have a radius of 0.5 pm, meaning that 1 LU= 0.4 pm (see above). The spacing between actin 
cables (or "lanes") was fixed to du = 2 pm= 5 LU. For computational reasons, it is desirable for the distance of 
particles from the wall to be such that there is at least one fluid node between wall and particle surface (otherwise 
lubrication forces need to be considered). Particle centres aie therefore placed at Zp = 2.35 LU (and they do not 
extend below z = 1.1 LU). The distance between the wall (mid-grid at = 0.5) and the particle surface is thus 
dw-p = 0.6 LU. This agrees well with the actual distance in the biological system: diameters of actin cables are 
between 100 nm and 200 nm (l^ and the size of Cham myosin is estimated to be about 100 nm 1 14], much larger 
than typical animal myosin motors. This results in about dw-p = 0.2 — 0.3 pm = 0.5 — 0.75 LU. Initial spacing 
between particles (centre to centre) varies with density: a separation d^.^ = 12.5 LU for instance corresponds to a 
line density of 0.2 (in units where 1 corresponds to a solid line of particles within a lane). 



We fixed the vacuole viscosity to r// = 0.02 and varied -qu in line with existing estimates; see Ref. Ill IP where 
the authors assumed a ratio rpi/i]! between 6 and 250. In simulations the vacuole has thickness 75LU=30//m. 
While this is much too small as the vacuole thickness can be ~ 500 ^m, the location of the upper boundary has 
virtually no effect on the velocities on the endoplasm-vacuole interface, which we are ultimately interested in. We 
also note that moving the upper vacuole boundary further away would, if anything, only increase the cytoplasm's 
velocity with respect to particle velocities. 

Finally, we note that here we assume that the viscosity of the thin layer in the cytoplasm which cannot be 
reached by the crowding agents is the same as that of the vacuole. This results in a slip layer of thickness equal 



"This is typical in LB work . Because the LB algorithm uses fluid inertia to update the dynamics, the most efficient discretisation 

uses the highest density compatible with still remaining in the low Reynolds number regime. With our chosen density the Reynolds number 
associated with the vesicle size never exceeds O.L 



to 220 nm, when p = 0.9. The thickness of the low-viscosity layer is not known but is restricted by the particle to 
wall distance dw-p = 0.6 LU, since the particles have to move within the high-viscosity region to provide effective 
drag on the main endoplasm layer. 



Results 

Fig. 1 shows a cartoon of the geometry of our simulations, which we recall consists in a two-layer fluid, modelling 
the endoplasm- vacuole system, with a slip length boundary condition, Eq. [6l A 2D carpet of molecular motors 
with their attached cargoes is arranged along regularly spaced actin cables and each of the motors is further subject 
to a constant external forcing which drives it along the lane. 



Allowing for a finite wall slip can reproduce the phenomenology of cytoplasmic streaming 




Figure 2: (Online version in colour.) (a) Ratio of average fluid velocity over average particle velocity u{z)/vp for 
a no-slip wall below the vesicles, (b) Partial slip at the bottom wall, coiTcsponding to a layer of lower viscosity and 
of thickness 220 nm (or 0.55 LU) close to the wall, can explain the high fluid velocities observed in experiments. 
Viscosity ratio is "quMi = 50, particle surface coverage = 0.1 with the distance between actin bundles being 
Ay = 5 LU corresponding to 2 /^m. The height of the cytoplasm is i^cyto = 25 LU=10 ^m. 

Fig. 2 shows the mean fluid velocity u{z) as a function of z, averaged over x and y and normalised by the 
average velocities of the motor carpet - the latter are recomputed for each simulation as they may vary. Fig. 2a 
considers no-slip boundary conditions, p = 0. Our direct simulations of the 3D flow within the endoplasm suggest 
that the velocities within the fluid sharply drop to about 50% of the motor-cargo velocities close to the wall. Inter- 
estingly, this is very much in line with the semianalytical calculations presented in Ref. [Tl] and relying on more 



assumptions than our LB simulations. As observed in Ref. Illll] . this is not compatible with experimental measure- 
ments which suggest that the velocity distribution within the endoplasm is within 10% of the motor velocities on 
the actin cables. The no-slip simulations also show a small peak in u{z), just above the particles. Increasing the 
size of the particles such a peak disappears, so this is probably a discretization artefact. 

Fig. 2b instead shows what happens if we allow for wall slip. As mentioned in the previous section, in the 
simplest approximation slip may be characterised by a single length scale, henceforth denoted as the slip layer 
thickness. This is equal to 220 nm in Fig. 2b, which is a reasonable value given the size of the cargoes, and 
corresponds to a slip parameter p = 0.9 in our LB simulations. The cytoplasmic flow is now dramatically different 



and much faster with respect to the no-shp case, and a tracer would move at over 90% of the mean motor speed, 
in agreement with the observations. 



(a) (b) 




Velocity w/vp Velocity u/vp 

Figure 3: (Online version in colour.) (a) Distribution of the fluid velocities inside the cytoplasm for no slip bound- 
ary conditions. The fluid velocities are normalised with respect to the average vesicle velocity. The distribution 
is sharply peaked, but at less than half the average vesicle velocity, (b) Same for a slip parameter p = 0.9. The 
distribution is now peaked at a value close to the average vesicle velocity in agreement with experiments. 

To further compare experiments and simulations, it is useful to go beyond averages and map out the velocity 
distribution in the cytoplasm. The distribution of velocities of organelles in the cytoplasm is available indirectly 
from light scattering experiments measuring Doppler shift in algal cells in vivo: these suggest that the distribution 
is sharply peaked with only about 10% variation close to the estimated myosin velocity Jsj]. Fig. 3 shows the 
distributions of the velocities in our simulations, computed on a regular array of points inside the endoplasm, 
which may be ascribed to tracer organelles within the cell. It can be seen that in the case where slip is allowed our 
simulated distribution is sharply peaked, and the standard deviation is about 10% or less of the mean value, in good 
semiquantitative agreement with experiments. In the no-slip case, the distribution of the endoplasmic velocities is 
still sharply peaked, but at a much lower velocity than that of the motors, at odds with experiments. 



The flow inside the endoplasm is quasi-one dimensional 

Our simulations consider a fully 3D flow, with a typical simulation volume equal to 25 x 20 x 100 in lattice units 
(see Choice of parameters section for a mapping to real units). As we directly model the two-dimensional caipet 
of moving motors, together with the interaction with the wall, it is reasonable to expect that close to the motors 
there will be significant inhomogeneities of the fluid velocity. This may be further increased by any hydrodynamic 
clustering of the motors. It is then useful to ask how deep inside the endoplasm the correlations and fluctuations 
in the fluid velocity imparted on the intracellular solvent by the motor dynamics persist, and what the decay of 
such correlations is. Fig. 4 shows a cut of the fluid velocity profile on different planes, stacked along the velocity 
gradient direction. 

As can be seen in Fig. 4, close to the particles the flow is indeed inhomogeneous, but surprisingly weakly 
so. The square root of the velocity variance is only about 1-10% of the typical fluid velocity even close to the 
motors; it then decays exponentially with distance and becomes negligible beyond 10 LU away from the walls. In 



other words, our simulations show that the flow is quasi-lD, con^oborating the analysis of Ref. llllll which assume 



this geometry at the outset, and also justifying the use of uniform velocity boundary conditions for the vacuole in 
Ref.H. 
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Figure 4: (Online version in colour.) (a) Decay of fluid velocity inhomogeneities with height z as relative variance 

o"u,rci = ulzf where i runs over all fluid nodes in a layer of fixed z. Inhomogeneities decay 

exponentially with z. Panels (b) to (d) show velocity maps at slices of fixed z for z = 2, z = 5 and z = 10 LU. 
Note the difference in scales for panels (b) to (d). 

Experiments constrain the model values of viscosities, motor density and slip layer thickness 

While Figs. 2 and 3 show that the key features of cytoplasmic streaming in plant cells may be reproduced by a 
reasonable set of values for cytoplasm viscosities, motor geometry, density and slip layer thickness, it should be 
noted that some of these control parameters are not known with high accuracy from experiments. It is therefore 
important to assess how our results and conclusions are affected if these ai^e changed, within a biologically relevant 
range. In this section we therefore analyse the impact of parameter variation on the fluid velocity within the 
endoplasm. 

We first consider the effect of the thickness of the slip layer, e. Fig. 2b referred to e = 220 nm, and Fig. 2a to 
e = 0, with only the former data accounting for a realistic streaming within the cytoplasm. Fig. 5 shows the effect 
of varying the slip layer, from 50 nm, which is smaller than the size of the motor, to 500 nm, the size of the cai^goes. 
Within this slip layer, we assume that the cytosol is basically an aqueous fluid - motivated by the view that the 
larger viscosity inside is essentially due to macromolecular crowding with crowding agents that do not reach up to 
the wall thus resulting in a lower viscosity close to the wall. The mean velocity shows a saturation behaviour as a 




Figure 5: (Online version in colour.) Influence of 
thickness of slip layer e with all other parameters 
{rjii/rji, (j), Ay and i/cyto) held constant. Thicknesses 
of over 200 nm for the low-viscosity layer are neces- 
sary to explain the high u{z)/vp ratio of over 90%. 



Figure 6: (Online version in colour.) Variation of sur- 
face coverage cj) and spacing of actin bundles Ay. Un- 
surprisingly, a denser coverage with particles entrains 
the cytoplasm to higher velocities. A surface coverage 
of 0. 1 is necessary to reach high enough velocities, all 
other parameters being equal. The distribution of parti- 
cles on the surface also plays a role: If the particles are 
dense in the x-direction but spaced further apart in the 
y-direction the ratio u{z)/vp is lower than if particles 
are denser in y and sparser in x (at constant overall </>). 
This effect is due to channels of slower fluid between 
actin bundles if they ai^e spaced further apart. 



function of e, with essentially any value larger than 200 nm enough to account for experimental values. The lower 
values for slip layer thickness, while geometrically consistent, lead to a velocity of the fluid which is about 80% 
of the motor velocity - this is slightly too low given current measurements, although the corresponding values of 
e cannot be definitely discarded. 

Fig. 6 shows the effect of changing the density of motors, measured by their surface coverage (p — ^occ / ^tot 
where ^occ is the surface area occupied by cargo beads and Atot is the maximum area available to them without 
blocking cai^goes on adjacent lanes (i.e. cai^goes placed on a rectangular lattice just touching their neighbouring 
cargoes on the same lane and on adjacent lanes). It is interesting to note that, once the slip layer thickness is large 
enough, even a surface density equal to 0.1, i.e. 10 times smaller than full coverage, is enough to account for the 
experimental velocities. Fig. 6 also shows that the distribution of particles on the surface also plays a role. For 
instance, one may achieve the same surface coverage = 0.1 by using fewer actin lanes with more myosin motors 
on each. The largest ratio between fluid and particle velocity is obtained when actin lanes are closely spaced. This 
effect is due to the creation of channels of slower fluid between actin lanes if the lanes are spaced further apart. 

Finally, Fig. 7 and Fig. 8 show the effects of cytoplasm viscosity and thickness respectively. The viscosity of 
the cytoplasm has been measured by several authors and values in the range 6-250 cP have been reported, whereas 
the viscosity of the vacuole may be taken to be close to that of water, 1 cP. The cytoplasm viscosity may actually 
depend on probe size and this makes it difficult to assess. For our case, as organelles were used in the original 
experiments to track the flow, it may be more appropriate to use viscosities observed with probes of abou the same 
size, hence towards the top end of the range quoted above. Fig. 7 shows that the cytoplasm viscosity does have a 
sizeable effect on the average cytoplasmic streaming velocity: the thr^eshold beyond which realistic streaming may 
be achieved via cargo dragging is about 20 cP, for which the average fluid velocity is about ~ 90% of the average 




Figure 7: (Online version in colour.) Variation of 
viscosity (ratio of viscosities). A viscosity ratio of 
7777/7/7=50 can account for fluid velocities above 90% 
of particle velocities while a ratio of 7/77/777=20 is 
not quite enough. For velocity ratios larger than 
7/77/777=100 u{z)/vp does not change much anymore. 
Note that for low viscosity ratios the velocity gradient 
in the cytoplasm also increases such that the velocity 
at the membrane becomes notably different from that 
at the bottom wall. 



Figure 8: (Online version in colour.) Plot of the mean 
velocity profiles as a function of height along z, for 
different thicknesses of the cytoplasm, -fTcyto- One may 
observe that the thickness of the cytoplasm plays no 
role for the ratio u{z)/vp as long as the top wall is 
in the far field and the viscosity of the vacuole much 
lower than that of the cytoplasm. Note that in the curve 
referring to //cyto = 20 /um the size of the vacuole has 
been increased to the minimal one which did not affect 
the fluid dynamics in the cytoplasm. 



motor velocity. Fig. 8 on the other hand shows that the thickness of the cytoplasm plays a comparatively much 
less important role. 



Discussion and conclusions 

The work we have presented provides a microscopic simulation of cytoplasmic streaming in a geometry relevant to 
that of plant cells, where there is an endoplasmic layer surrounding a much larger vacuole, and the flow inside the 
endoplasm is driven by molecular motors which run on actin lanes dragging vesicles. Therefore our calculations 
can probe the viability of current continuum fluid dynamic theories of cytoplasmic streaming, which often lump the 
microscopic detail of the motor motion into a shear velocity boundary condition for a purely continuum problem 
inside the vacuole Isj]. 

Our main conclusion is that if motors on the cytoplasm are attached to vesicles or organelles, the drag exerted 
on the fluid through their motion is enough to lead to an essentially continuum flow within the endoplasm and the 
vacuole, with average fluid velocity fully consistent with those observed in cytoplasmic streaming. This conclusion 
holds for a range of realistic pai^ameters, and suggests that there is no need to invoke a solid-like coupling between 
the motors moving along the subcortical actin tracks and the cytoskeletal network, or endoplasmic reticulum. 
Such an elastic coupling was deemed necessary in order to account for the experimental velocities in the previous 



work of Ref. lilll] . and it is therefore useful to point out some key differences between our approaches. Most 
important is the fact that our work presents direct hydrodynamic simulations of the motion of solid spheres close 
to a boundary in a two-layered fluid mimicking the endoplasm- vacuole of an algal cell. As such, we fully take into 
account sphere-sphere and sphere-wall interactions, and crucially also consider slip at the wall. The hydrodynamic 
treatment in Ref. Illltl on the other hand neglected these interactions, and this is probably at the origin of the 



quantitative discrepancy between our conclusions (although it should be noted that even the results of Ref. Ill ill 



did not definitely rule out the purely viscous coupling advocated in our work). Indeed, many-body hydrodynamics 
leads to a collective effect which increases the drag exerted on the fluid by a moving carpet of motors, whereas the 
wall, which could counteract that effect in principle, does not counteract it in our simulations due to the inclusion 
of a slip boundary condition. Our approach of introducing slip boundary conditions can likely be generalised 
to a range of biological solvents whose high effective viscosity is due to macromolecular crowding, so that a 
layer of lower effective viscosity forms close to walls, where the crowding agents (which are macromolecules and 



It is also important to recognise some remaining limitations of our calculation. Firstly, the interior of the 
plant cell we wish to model is actually cylindrical, while we have modelled a slab of fluid. This is perhaps not 
too important as the radius of curvature of the cell is much larger than the thickness of the endoplasm, which is 
the region of interest, as this is where the coupling between motor motion and vacuolar fluid dynamics occurs. 
Secondly, we have neglected all details of the intracellular fluid and just modelled the cytoplasm/endoplasm as a 
Newtonian viscous fluid. In reality, the cytoplasm is non-Newtonian, shear thinning, and it would be interesting to 
see how a more careful model of the intracellular solvent may affect our results. To this end, we should however 
make a decision as to whether to model the endoplasm as simply a crowded medium or to include the endoplasmic 
reticulum as a polymer network. Given the uncertainty with which several of the pai^ameters we use are known, 
this does not yet seem appropriate and might have to await more detailed experimental study of the rheological 
propertites of the endoplasm. In this context it is also interesting to note that in Ref. where the authors studied 
cytoplasmic streaming in C. elegans, simulations at constant viscosity proved sufficient to quantitatively reproduce 
the streaming motion of the (actually non-Newtonian) cytoplasm. A final simplification of our treatment is that 
the membrane separating endoplasm and vacuole was computationally treated as an infinitesimally thin interface 
where the fluid velocity was continuous, whereas in reality this membrane is better described as a viscoelastic 
membrane. 

This work was supported by a DAAD Postdoctoral Fellowship to KW and by EPSRC Grant EP/E030173. 
MFC holds a Royal Society Professorship. 

References 

[1] Squires, T. M. 2010. A furtive stare at an intra-cellular flow. /. Fluid Mech. 642, 1^. 

[2] Shimmen, T. 2007. The sliding theory of cytoplasmic streaming: fifty years of progress. /. Plant. Res. 120, 



[3] Shimmen, T. and Yokota, F. 2004. Cytoplasmic streaming in plants. Curr. Opin. Cell Biol. 16, 68-72. 

[4] Niwayama, R., Shinohara, J., and Kimura, A. 2011. Hydrodynamic property of the cytoplasm is sufficient 
to mediate cytoplasmic streaming in the caenorhabiditis elegans embryo. Proc. Natl. Acad. Sci. USA 108, 



[5] van de Meent, J.-W., Tuval, I., and Goldstein, R. F 2008. Nature's microfluidic transporter: Rotational 
cytoplasmic streaming at high Peclet numbers. Phys. Rev. Lett. 101, 178102^. 

[6] Goldstein, R. F., Tuval, I., and van de Meent, J.-W. 2008. Microfluidics of cytoplasmic streaming and its 
implications for intracellular transport. Proc. Natl. Acad. Sci. USA 105, 3663-3667. 

[7] van de Meent, J.-W., Sederman, A. J., Gladden, L. F., and Goldstein, R. F. 2010. Measurement of cytoplasmic 
streaming in single plant cells by magnetic resonance velocimetry. /. Fluid Mech 642, 5-14. 




31-43. 



11900-11905. 



[8] 



Mustacich, R. V. and Ware, B. R. 1974. Observation of protoplasmic streaming by laser-light scattering. 
Phys. Rev. Lett. 33, 617-620. 



[9] Mustacich, R. V. and Ware, B. R. 1976. A study of protoplasmic streaming in nitella by laser doppler 
spectroscopy. Biophys. J. 16, 373-388. 

[10] Sattelle, D. B. and Buchan, R 1976. Cytoplasmic streaming in chara coralUna studied by laser light scattering. 
J. Cell Sci. 22, 633-643. 

[11] Nothnagel, E. A., and Webb, W. W. 1982. Hydrodynamic models of viscous coupling between motile myosin 
and endoplasm in characean algae. /. Cell. Biol. 94, 444-454. 

[12] Kachar, B., and Reese, T. S. 1988. The mechanism of cytoplasmic streaming in characean algal cells: sliding 
of endoplasmic reticulum along actin filaments. J. Cell Biol. 106, 1545-1552. 

[13] Yamamoto, K., Shimada, K., Ito, K., Hamada, S., Ishijima, A., Tsuchiya, T., and Tazaw, M. 2006. Chara 
myosin and the energy of cytoplasmic streaming. Plant Cell Physiol. 47, 1427-1431. 

[14] Morimatsu, M., Nakamura, A., Sumiyoshi, H., Sakaba, N., Taniguchi, H., Kohama, K., and Higashi-Fujime, 
S. 2000. The molecular structure of the fastest myosin from green algae, Chara. Biochem. Biophys. Res. 
Commun. 270, 147-152. 

[15] Cates, M. E., Stratford, K., Adhikari, R., Stansell, R, Desplat, J. C, Pagonabarraga, I., and Wagner, A. J. 
2004. Simulating colloid hydrodynamics with lattice boltzmann methods. /. Phys. Condens. Matt. 16, 
S3903-S3915. 

[16] Nguyen, N. Q. and Ladd, A. J. C. 2002. Lubrication corrections for lattice-boltzmann simulations of particle 
suspensions. Phys. Rev. E 66, 046708. 

[17] Sbragaglia, M. and Succi, S. 2005. Analytical calculation of slip flow in lattice boltzmann models with 
kinetic boundary conditions. Phys. Fluids 17, 093602. 

[18] Henrich, O., Stratford, K., Marenduzzo, D. and Cates, M. E. 2010. Ordering dynamics of blue phases entails 
kinetic stabilization of amorphous networks. Proc. Natl. Acad. Sci. USA 107, 13212-13215. 

[19] Nagai, R. and Rebhun, L. I. 1966. Cytoplasmic microfilaments in streaming nitella cells. /. Ultrastruct Res. 
14, 571-89. 

[20] Eisenriegler, E. 1993. Polymers Near Surfaces. World Scientific, Singapore. 



Appendix 

In this Appendix we give additional information on our numerical framework and on the derivation of the formulas 
determining the partial shp parameter in our simulations, see Eq. [5l|7J 

In order to model slip at a solid wall in lattice Boltzmann (LB) simulations, a partial slip parameter p can 
be introduced which specifies the reflectivity of the wall. This parameter determines the fraction of a velocity 
population which is reflected specularly (or bounced forward) at the wall in an LB streaming step (see Fig. |9ll. 
Note that p = corresponds to full no slip boundary conditions for the velocity field at the wall, whereas p = 1 
corresponds to full slip - we refer to any value of p in between as partial slip. The wall reflectivity can be 



interpreted physically in terms of the Knudsen number which enters the simulation via the relaxation time r Ill7ll . 
Another possibility is to interpret the wall slip in a more general way where p is simply an open pai^ameter that 
can be determined if the slip length n(0)/n'(0) is known from experiment or theory. 

One particular setting that can be described in this way is that of a two-viscosity system where a thin low- 
viscosity layer is situated close to the wall and a bulk fluid of higher viscosity beyond. If the behaviour of the 



before streaming 



after streaming 



bounce back bounce forward 

fraction q = l - p fraction p 




Figure 9: (Online version in colour.) Velocity populations impinging on the wall are partially reflected forwards 
(fraction p) and partially bounced back (fraction q = 1 — p). 



thin layer itself is not of interest, all the layer does is introduce effective slip for the bulk system which can be 
simulated using the slip parameter p. 

Assuming the velocity profile of the fluid to be of the form 



E^=o 0<z<e 



u{z) 



the slip length u{0)/u'{0) can be determined analytically for e <^ R where R is the total system size and e the 
thickness of the low-viscosity layer. Various conditions on u give 

u(0) = =^ Oo = no slip at wall 

N 

ii(e+) = u(e_) =^ bo = ''^{an — bn) = {ai — bi)e + 0{e^) continuity of m 

n.=l 

rjiiu'U^) = Tjiu'U^) ai = — bi + Oie) continuity of stress 

and therefore 



bo = b,[J^-l]e + 0{e^) 



Observing that u(0) = 6o and li'(O) = 6i we arrive at 



1 (8) 



n'(o) V m 

to first order in e, where r]i is the (lower) viscosity of the thin layer close to the wall and r?// the (higher) viscosity 
of the bulk system. 

In order to relate slip length n(0)/n'(0) and slip parameter p and relaxation time r in the LB method, we 
consider the D2Q9 model. All simulations have been performed in the D3Q19 model where the derivation works 
in exactly the same way but the algebra is more cumbersome (Fig.fTOl). 

In a collision step, velocity populations of the same node collide and relax towards their local and instantaneous 
equihbrium population f1'^{t, x) 

/.(f,x) = /*(t,x)-i(/*(i,x)-/ni,x)) 

= /*(t,x)('l-i]+i/r(t,x) (9) 



15 




Figure 10: (Online version in colour.) Left: Set of two-dimensional velocities ej with i = to 8, D2Q9, used 
in the following derivations. Right: Set of three-dimensional velocities with i = to 18, D3Q19, used in the 
simulations. 

Here fi {t, x) is the velocity population of direction i at time t and node x after collision and /* (t, x) the population 
before collision. The relaxation time r determines how far the population relaxes in each step and is related to 
viscosity by r = . Here only the bulk viscosity r/// enters the formula as it is the bulk fluid that is simulated 

in LB whereas the low-viscosity layer and thus r// is incorporated into the boundary condition. 
The equilibrium populations are given by 

fPiz) = pcji (^1 + 3ei ■ u(t, x) + ^ (e, • u(t, x))^ - ^u(t, x)^^ (10) 

with oji = 0J3 = 0J5 = LJj = 1/9 =: wi and 002 = u}4 = ujq = ujg = 1/36 =: W2- Density p is assumed to be 
constant and fluid velocity is given by 

pu{t,x) = ^eifi{t,x). (11) 

i 

We next need to make a couple of assumptions on the velocity profile. First, we assume that u(t, x) = 
u{t, z)ex is parallel to the wall and that it only depends on the distance to the wall z. So for the velocity in 
x-direction we have 

pu{z) = fl{z) + fi{z) - fliz) - ft{z) - n{z) + f;{z) (12) 

In the streaming step velocity populations from time step t are passed on to neighbouring nodes and become 
pre-collision populations of time step t + 1. Below, streaming is given for the layer of fluid nodes just above the 
wall, i.e. those experiencing reflection from the wall. As velocities are assumed to be homogeneous in x- and 
y-directions those indices are omitted altogether and /j ^ denotes velocity population i in layer k above the first 



fluid layer. For the first layer k 



will be omitted, 

flit + I 

flit + 1 
f^it+i 



flit) 

qfeit) + Pfsit) 
frit) 

Pfeit) + qfsit) 

/5(i) 

feAt) 

frAt) 

fs,iit) 



(13) 



The next assumption to be made is that steady-state has been reached and the velocity profile and all populations 
do no longer change with time meaning that all reference to t can be dropped. 
Plugging streaming (Eg. [T3]) and collision (Eq.|9l) into Eq.[T2]yields 
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To close this equation we need to express the pre-coUision populations /* in terms of equilibrium populations f^'^ 
as for those we have Eq. [10] relating them to p and u. 

For populations /f and /g this is easy as no populations from higher layers (larger z) are involved and the 
profile is assumed to be homogeneous in x meaning that after many iterations fl and /g simply relax to f^^ and 
/g^. For /| and /| this is more complicated as populations from higher and higher layers wiU be involved 



f*,it + n) = flnit)[l 



^ k=l 
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6,k- 



(14) 



For large enough n the term involving 
the equation for pu leaves us with 



PUQ 



1 — i)" can be disregarded and using the above expression and Eq.fTOlin 
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where uq is the velocity in the layer of fluid nodes closest to the wall and the velocities in higher layers. 

Using Taylor expansions Uk = 
Uk and Uk+i in the sum above finally results in 



°° 1 jtUq^ for velocities at higher layers and combining subsequent terms of 



Uq = Wi6uq H W2G 
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(15) 



In principle ^^^=0 ~ t)^ evaluated for arbitrary j but then more and more derivatives Uq^ have to be 

known in advance. We therefore truncated after the linear term and checked that the resulting formula still holds 
within our numerical precision for Poiseuille flow where the velocity profile in the channel is quadratic. 
In the linear approximation Uq =0 for j > 1 we obtain 

Uq = wiGuo + -W26 {t2puq + (r^ - r) 2pu'q + tu'q) 

2 1 1,1, , , 

= 3^^o + -^puo + g^O + 3 ("^ - 1) P^Q- 

At this point it is important to note that in line with standard LB algorithms working within the Euler scheme which 
we use, the wall is implemented between two nodes, so that it is mid-grid. Therefore ti(0) from Eq. [5]does not 
con^espond to the velocity at the first layer of fluid nodes uq but lies half a lattice spacing below n(0) = uq — 1/2u'q. 
As Uq is assumed constant it is of course the same as u'(0) and we get 

u{0) _ 2t - 1 p 

u'{0) 2 1-p ^ ' 



which can be used together with Eq.[8]to obtain the formula relating slip parameter p and physical quantities e, r]i 
and r/// in Eq.|7]of the Methods and models section. 



